﻿#include <amp.h>
#include <amp_math.h>

using namespace concurrency;
using namespace concurrency::fast_math;

#define ⇶(name) namespace name {
#define ⟡ }
#define → restrict(amp, cpu) ->
#define def auto
#define ← =
#define let auto
#define ⩾ >=
#define ≈ ==
#define √ sqrtf
#define ⦁ ,
#define import(name) using namespace name;
#define ∙ *
#define ≺
#define ≻
#define ∷ =
#define ≔ { return
#define enddef };

typedef float ℝ;

⇶(structs)
template <class T, int dims>
struct Κ {
	T elements[dims];
	enum xxx {dimensions = dims};

	let operator[] (int index) → T&	{
		return elements[index];
	}

	let operator[] (int index) const → T const& {
		return elements[index];
	}

	let operator+= (Κ<T, dims> const& other) → Κ<T, dims>& {
		for (unsigned int index = 0; index < dims; ++index) {
			(*this)[index] += other[index];
		}

		return *this;
	}
};

typedef Κ<float, 2> ℝ²;

template <class U, class T>
struct Φ {
	Κ<T, U::dimensions> position;
	T height;
	T width;
};
}

⇶(functions)
	import(structs)

	struct pseudogauss {
	/*                                                                                                                */
	def static f(Φ<ℝ², ℝ> const& potential ⦁ ℝ² const& κ) → ℝ {
		let Δx ← (κ[0] - potential.position[0]);
		let Δy ← (κ[1] - potential.position[1]);
		let ‖Δ‖² ← Δx * Δx + Δy * Δy;

		if (‖Δ‖² ≈ ℝ(0))
			return ℝ(potential.height);

		let width² ← potential.width * potential.width;

		if (‖Δ‖² / width² ⩾ ℝ(1))
			return ℝ(0);

		let ‖Δ‖∕width ← √(‖Δ‖²) / potential.width;

		let ⦗‖Δ‖∕width⦘² ← ‖Δ‖∕width * ‖Δ‖∕width;
		let ⦗‖Δ‖∕width⦘³ ← ‖Δ‖∕width * ⦗‖Δ‖∕width⦘²;

		return (ℝ(2) * ⦗‖Δ‖∕width⦘³ - ℝ(3) * ⦗‖Δ‖∕width⦘² + ℝ(1)) * potential.height;
	enddef
	/*                                                                                                                    */
	def static ∂f∕∂Κ_x(Φ<ℝ², ℝ> const& potential ⦁ ℝ² const& κ) → ℝ {
		let Δx ← (κ[0] - potential.position[0]);
		let Δy ← (κ[1] - potential.position[1]);
		let ‖Δ‖² ← Δx * Δx + Δy * Δy;

		if (‖Δ‖² ≈ ℝ(0))
			return potential.height;

		let width² ← potential.width * potential.width;

		if (‖Δ‖² / width² ⩾ ℝ(1))
			return ℝ(0);

		let ‖Δ‖∕width ← √(‖Δ‖²) / potential.width;

		return Δx * ℝ(6) * (‖Δ‖∕width - ℝ(1)) * potential.height;
	enddef
	/*                                                                                                            */
	def static ∇_Κf(Φ<ℝ², ℝ> const& potential ⦁ ℝ² const& κ) → ℝ² {
		ℝ² c;

		let Δx ← (κ[0] - potential.position[0]);
		let Δy ← (κ[1] - potential.position[1]);

		float d = Δx * Δx + Δy * Δy;

		if(d == 0.0f)
		{
			c[0] = 0.0f;
			c[1] = 0.0f;
		}
		else
		{
			d = sqrtf(d) / potential.width;

			if( d>=1.0f )
			{
				c[0] = 0.0f;
				c[1] = 0.0f;

				return c;
			}

			float p = 6.0f*(d-1.0f)*potential.height/potential.width;
			c[0] = Δx*p;
			c[1] = Δy*p;
		}

		return c;
	enddef
};
	template <class function>
	struct ∑
	{
		def static f(array<Φ<ℝ², ℝ>, 1>& Φs, ℝ²& κ) → ℝ
		{
			auto sum = ℝ(0);

			for (unsigned int index = 0; index < Φs.extent.size(); ++index)
				sum += function::f(Φs[index], κ);

			return sum;
		}

		def static ∇_Κf(array<Φ<ℝ², ℝ>, 1>& Φs, ℝ²& κ) → ℝ²
		{
			auto sum = ℝ²();

			for (unsigned int index = 0; index < Φs.extent.size(); ++index)
				sum += function::∇_Κf(Φs[index], κ);

			return sum;
		}
	};
}